###Generic Function for T Rejection Point
T.Rejection.Point <- function(alpha,df,tails){
  if(tails==2)return(qt(1-alpha/2,df))
  if((tails^2) != 1) return(NA)
  return(tails*qt(1-alpha,df))
}
### Generic Function for T-Based Power
Power.T <- function(delta,df,alpha,tails){
  pow <- NA
  R <- T.Rejection.Point(alpha,df,abs(tails))
  if(tails==1)
    pow <- 1 - pt(R,df,delta)
  else if (tails==-1)
    pow <- pt(R,df,delta)
  else if (tails==2)
    pow <-  pt(-R,df,delta) + 1-pt(R,df,delta)
  return(pow)
}
### Power Calc for One-Sample T
Power.T1 <- function(mu,mu0,sigma,n,alpha,tails){
  delta = sqrt(n)*(mu-mu0)/sigma
  return(Power.T(delta,n-1,alpha,tails))
}


### Plot Power Curve
curve(Power.T1(75,70,10,x,0.05,1),
      10,100,xlab="Sample Size",
      ylab="Power",col="red")

#### Home in
curve(Power.T1(0.80,0,1,x,0.01,2),
      10,100,xlab="Sample Size",
      ylab="Power",col="red")
abline(h=.95)
abline(v=30)
abline(v=35)

Power.T1(0.80,0,1,32,0.01,2)

curve(Power.T1(0.80,0,1,x,0.01,2),
      30,35,xlab="Sample Size",
      ylab="Power",col="red")
abline(h=.95)
abline(v=31)
abline(v=32)

Power.T2 <- function(mu1,mu2,sigma,n1,n2,alpha,
                     tails,hypo.diff=0){
  delta = sqrt((n1*n2)/(n1+n2))*
    (mu1-mu2-hypo.diff)/sigma
  return(Power.T(delta,n1+n2-2,alpha,tails))
}

Power.T2(0.50,0,1,20,20,0.05,2)



Power.GT <- function(mus,ns,wts,sigma,alpha,
                      tails,kappa0=0){
     W = sum(wts^2/ns)
     kappa = sum(wts*mus)
     delta = sqrt(1/W) * (kappa-kappa0)/sigma
     df = sum(ns)-length(ns)
     return(Power.T(delta,df,alpha,tails))
 }

Power.GT(c(75,75,70),c(10,10,10),c(1,1,-2),10,0.05,2)

Power.T2Correlated <- function(mu1,mu2,sigma,n,rho,alpha,tails){
  delta = sqrt(n)*(mu1-mu2)/sigma/sqrt(2*(1-rho))
  return(Power.T(delta,n-1,alpha,tails))
}
Power.T2Correlated(15,10,10,100,0.50,0.05,2)


curve.js <- function(f,a,b,points=100,type='l',...){
  ftext <- paste("g <- function(x){",f,"}")
  eval(parse(text=ftext))
  x <- seq(a,b,length=points)
  plot(x,mapply(g,x),type,...)
}
plot.curve <- function(f,a,b,points=100,type='l',...){
  x <- seq(a,b,length=points)
  plot(x,mapply(f,x),type,...)
}